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\^ ' Abstract. Recently we have found that a family of models of partially relaxed, anisotropic stellar systems, inspired earlier 

, by studies of incomplete violent relaxation, exhibits some interesting thermodynamic properties. Here we present a systematic 

T-H ■ investigation of its dynamical characteristics, in order to establish the basis for a detailed comparison with simulations of 

^ ' collisionless collapse, planned for a separate paper. For a full comparison with the observations of elliptical galaxies, the 

0^ , models should be extended to allow for the presence a sizable dark halo and of significant rotation. In the spherical limit, the 

' family is characterized by two dimensionless parameters, i.e. f, measuring the depth of the galaxy potential, and v, defining 

the form of a third global quantity Q, which is argued to be approximately conserved during collisionless collapse (in addition 
to the total energy and the total number of stars). 

The family of models is found to have the following properties. The intrinsic density profile beyond the half-mass radius fm 
, is basically universal and independent of 4'. The projected density profiles are well fitted by the law, with n ranging from 

■ 2.5 to 8.5, dependent on with n close to 4 for concentrated models. All models exhibit radial anisotropy in the pressure 

' tensor, especially in their outer parts, already significant at r x r^. At fixed values of v, models with lower T are more 

anisotropic; at fixed values off, models with lower v are more concentrated and more anisotropic. When the global amount of 
^ ■ anisotropy, measured by 2K, /Ky, is large, the models are unstable with respect to the radial-orbit instability; still, a wide region 

^ ' of parameter space (i.e., sufficiently high values of for v > 3/8) is covered by models that are dynamically stable; for these, 

, the line profiles (line-of-sight velocity distribution) are Gaussian at the 5% level, with a general trend of positive values of /?4 at 

radii larger than the effective radius R^. 
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1 . Introduction simulations of collisionless collapse) had been taken in terms of 

the so-called /„ models (Bertin & Stiavelli 1984 1. These mod- 

A simple picture for the incomplete violent relaxation of stellar constructed from a distribution function which, in the 

systems considers the collapse of a dynaimcally cold cloud of spherical limit, reduces to U = A(-Ef'^ exp (-aE - cfi 12), 

stars or star-clumps initially far from equilibrium. The study of ^^^^^^^^ ^^j^^^ v^nx^h^^ for positive values of E; 

this Ideal and relatively simple p rocess, pioneered by important ^^^^ ^ p^^jji^^ constants and E and J denote the 

analysis in the 60s (starting with|Lynden-Bell 1967t and simu- ^^^^^^^ ^^^^ ^^^^^^ ^„g^l^ momentum, 
lations in the 80s (see 'van Alba da 1982> . is still incomplete. A 

proper understanding of such process is a prerequisite to more From those earlier investigations, it immediately became 

ambitious attempts at constructing physically justified models clear that the physical clues gathered from the picture of coUi- 

of elliptical galaxies in which the problem of galaxy formation sionless collapse (i.e., that incomplete violent relaxation leads 

is set in the generally accepted cosmological context of hierar- to systems that ai-e well relaxed in their inner regions, r <sc r^, 

chical clustering. and characterized by radially anisotropic envelopes for r » km) 

Before the development of such cosmological scenarios, a do not lead uniquely to the models, but instead identify 

first step (in the direction of incorporating in a simple analyti- a wide class of atti-active distribution functions of which 

cal and physically justified framework the clues gathered from represents just one simple and interesting case. At that time, 

the primary goal of a series of investigations (described, e.g., 

Send offprint requests to: m.trenti@sns.it bv .Berlin & Stiavelli 1993 J was to test whether the models in- 
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spired by studies of collisionless collapse were realistic and 
thus could serve as a useful tool to interpret the observations. 
Indeed, the foa models turned out to "explain" the lumi- 
nosity law dde Vaucouleurs 1948 l and, extended to the two- 
component case, were used successfully to probe the presence 
and size of dark halos in elliptical galaxies. 

Because of the focus on such astronomical applications, 
the problem of a detailed comparison between the or other 
models and the products of collisionless collapse from N-body 
simulations was given lower priority and basically left aside. 
Yet, it was pointed out that, with respect to the products of 
numerical simulations of collisionless collapse, the /co models 
had the undesired feature of being too isotropic. Some authors 
( Merritt, Tremaine, & Johnstone 1989t thus suggested that the 
/co models should be extended to and used in the parameter 
domain where a < 0; however, this attempt failed, not only 
because a proper physical justification was lacking, but espe- 
cially because such "negative-temperature" models suffer from 
the opposite difficulty, i.e. they were shown to be so anisotropic 
that they are violently unstable and bound to evolve on a time 
scale even smaller than the typical crossing time. 

In a recent paper (IBertin & Trenti 2003J we have re- 
visited the problem of the structure and dynamics of par- 
tially relaxed stellar systems starting from a thermodynamic 
description. In view of the paradigm of the gravothennal 
catastrophe (see lAntonov 19621 |Lynden-BeIl & Wood 1968| 
IKatz 1978 1979 19801, we found it appropriate to approach 
the problem in terms of the so-called f^^^ models (for a defi- 
nition, see following Sect. 2); with respect to the fco models, 
the /^''^ models follow from a statistical mechanical deriva- 
tion JStiavelh & Bertin 1987 1 that is formally straightforward, 
while they are also known to have generally similar and reason- 
ably realistic structural properties. It should be noted that the 
dynamical properties of the /^''' models had not been studied 
earUer in detail. Such recent investigation, focusing on the is- 
sue of the gravothermal catastrophe, convinced us that the f^'^'' 
models do have many attractive features; in particular, these 
models turn out to have a higher degree of radial pressure 
anisotropy with respect to the foe models, making them more 
suitable to describe the results of simulations of collisionless 
collapse. Such encouraging preliminary inspection was the ba- 
sis for a thorough study that we have thus performed and that 
we present here. 

In this paper we provide a systematic description of the dy- 
namical properties of the models. In the spherical limit, 
these models define a two-parameter family, characterized by 
the dimensionless parameter ^F, measuring the depth of the 
galaxy potential, and the dimensionless parameter v, defining 
the form of a third global quantity Q (which is argued to be 
approximately conserved during collisionless collapse, in ad- 
dition to the total energy and the total number of stars; see 
Sect. 2 for the relevant definitions). In Sect. 3 we present in- 
trinsic and projected density profiles and fit the latter structural 
characteristics in terms of the R^^" law (Sersic 19681. We then 
proceed to illustrate in detail their phase space properties, in 
particular by calculating the relevant pressure tensor profiles, 
the pressure anisotropy content, the projected velocity disper- 
sion profiles, the line profiles (line-of-sight velocity distribu- 



tion), and the relevant phase space densities N{E, fi) and A^(£') 
(Sect. 4). In Sect. 5 we examine the stability of the models with 
respect to the radial orbit instability, not only by inspection of 
well-known global anisotropy indicators, but also with the aid 
of a number of N-body simulations. In Sect. |6lwe perform a 
first test of the viability of these models by applying them to 
the observed properties of NGC 3379, a galaxy that is charac- 
terized by lack of significant rotation and is likely to possess 
only small amounts of dark matter. The conclusions are drawn 
in Sect. 7. 

This paper also sets the basis for a detailed comparison 
between the /'''^ family of models and the products of colli- 
sionless collapse resulting from N-body simulations, to be pre- 
sented in a follow-up paper For a full comparison with the ob- 
servations of elliptical galaxies, the models should be extended 
to allow for the presence of a sizable dark halo and of signifi- 
cant rotation. 

1.1. Models of galactic structure: physical outlook 

The need for anisotropic models to describe elliptical galax- 
ies can be traced back to the strong empirical arguments that 
come from the observations of non-spherical geometry in the 
absence of significant rotation (e.g., see the discussion pro- 
vided by Bertin & Stiavelli 1993 and references therein); to be 
sure, direct evidence for pressure anisotropy in round systems 
is not easy to obtain, because so far the observed line pro- 
files show only modest deviations from a Gaussian (e.g., see 
IGerhard et al. 2001> . 

The main goal of this paper is to study the properties of a 
given specific dynamical context able to provide physical jus- 
tification for the existence of anisotropic equilibria, that is the 
dynamical framework of collisionless collapse. The process of 
collisionless collapse, together with its accompanying mecha- 
nism of incomplete violent relaxation, is one (but not the only) 
element that is expected to play an important role in the for- 
mation of stellar systems. In spite of the many papers that have 
addressed issues related to such dynamical context (in addi- 
tion to the papers cited earlier in the Introduction, e.g. see Shu 
1978, 1987, VogUs 1994, Hjorth & Madsen 1995), it is not yet 
clear whether an analytically tractable distribution function, or 
family of distribution functions, can be assigned to the prod- 
ucts of such collisionless collapse. This paper tries to provide a 
contribution to this problem. At the same time, we still need to 
establish how far the properties of such products, often studied 
in a simplified one-component picture, would be from those of 
observed stellar systems, and in which way. 

An enormous amount of work has focused and is currently 
focusing on the demands (on galactic structure) from the cos- 
mological context. In particular, this has led astronomers to 
look for the presence, in observed objects, of cuspy density 
distributions of dark matter, following a universal profile sug- 
gested by cosmological simulations (see Navarro, Frenk, & 
White 1997; Moore et al. 1998, Ghigna et al. 2000). Other 
studies have addressed the issue of the establishment of the 
galaxy scaling laws (such as the Fundamental Plane for ellip- 
tical galaxies; e.g., see Meza et al. 2003, Lanzoni et al. 2004) 
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in the accepted cosmological scenario. These demands are be- 
yond the scope of the present paper, but should eventually be 
faced, as a point of contact between dynamical investigations 
of individual galaxies and studies of the evolving universe from 
which they were formed. 

In this respect, a point to be noted, to avoid unnecessary 
confusion about the aims of purely dynamical studies, is the 
following. A priori, studies of collisionless collapse within the 
line of research adopted in this paper have nothing to say 
about some important issues such as the establishment of the 
Fundamental Plane (but see Gonzalez-Garcia & van Albada 
2003, Nipoti et al. 2003), because the relevant scaling laws de- 
pend on physics that goes beyond pure dynamics, which is in- 
herently scale-free. In turn, pure dynamical studies can try to 
explain why galaxies prefer the R^^* law, which is a structural 
property, instead of just accepting it as an empirical fact (as 
often done in a number of otherwise important astrophysical 
studies). 

Thus we would like to emphasize that this paper represents 
only one step in the direction of a comparison with the obser- 
vations. To deal fruitfully with the presence of dark matter and 
other important ingredients (such as significant rotation and the 
possible presence of an additional disk component), one first 
has to master the properties of one-component models, which, 
as shown in this paper, turn out to exhibit a variety of interest- 
ing dynamical properties. In fact, it is rewarding and a priori 
unexpected to find that, as a result of a simple conjecture about 
the way to characterize incomplete violent relaxation (the ad- 
dition of Q to the natural constraints under which Boltzmann 
entropy is extremized; see Sect. 2), one-component spherical 
models are identified able to fit products of N-body simulations 
over nine orders of magnitude in density (see Bertin & Trenti 
2004) and, at the same time, the observed photometric profile 
(over about ten magnitudes) and the inner kinematic profile 
(inside Re) of the best studied elliptical galaxy (see Sect. 6). 
Therefore, in spite of its incompleteness, the stage reached so 
far is definitely interesting from the physical point of view. 
Discrepancies with respect to the observations play the wel- 
come role of providing concrete indications about the role of 
the ingredients that are a priori ignored by the purely dynami- 
cal and highly simplified picture considered in this paper. 

2. Model construction and the relevant parameter 
space 

In the spherically symmetric limit, in order to allow for the pos- 
sibility that a stellar system is only partially relaxed, one may 
extremize the Boltzmann entropy S — - J fin fd^xd^w under 



the constraint that the total energy E/oi — (1/3)/ Efd^xd" 
the total mass M - ^ fd^xd^w, and the additional quantity 



w. 



Q 



\E\-^''''^fd\d^w 



(1) 



are constant jStiavelli & Bertin 1987t . Here E - /2 + <i> and 
7^ = |r X vfp represent the specific energy and the specific an- 
gular momentum square of a single star subject to a spherically 
symmetric mean potential <i>{r). (We have decided to use the 



symbol w, instead of v, for the velocity variable so as to avoid 
confusion with the symbol v.) If only E,ot and M were kept 
fixed, the extremization process would lead to a Maxwellian, 
that is to an isothermal and isotropic distribution function, ap- 
propriate for a fully relaxed system. Some arguments have been 
provided as to why a quantity such as Q should be, at least 
approximately, conserved iStiavelh & Bertin 1987 ) and those 
will not be repeated here. In any case, the conservation of Q 
should be taken as a conjecture. In a follow-up paper (Trenti, 
Bertin, & van Albada in preparation) we will discuss this prob- 
lem further by examining a set of collapse simulations, for 
which the conservation of Q will be tested. 

Such extremization leads to the following expression for 
the distribution function: 



= A exp 



-aE — d 



v/2 



(2) 



where v, a. A, and d are positive real constants. This set of 
constants provides two dimensional scales (for example a mass 
scale Mscaie = Afl^^^'^c/"-'^'' and a reference radius Rscaie - 
fl^'^^'c/"'^'') and two dimensionless parameters. For the latter 
two quantities, we may refer to v and y - ad^^^ liAnGA). The 
distribution function is taken to vanish for unbound orbits, that 
is for £■ > 0. 

The two-parameter family of models is then constructed by 
solving the Poisson equation: 



V2(l)(r) = 47rG ^ f'\r,w)d^w. 



(3) 



for the potential <l>(r). As briefly described by Bertin & Trenti 
(2003), the solution is obtained numerically after the equation 
has been cast in dimensionless form: 



1 d 



1 



r^ dr dr 



y 



(4) 



with f - r I Rscaie and <!) - a<b (for the adopted scaling, see the 
Appendix). The Poisson equation is integrated with "initial" 
conditions d^/dr = and <1) = at r = and the parameter 
y is determined as an eigenvalue, y = yi^), in order to satisfy 

the condition of Keplerian decay (<1) 1 /r) at large radii. The 

general behavior of the function yCi') is similar to that of the 
corresponding function for the foa models; after some oscilla- 
tions, 7 tends to a "plateau" at large values of *P (as indicated 
by Fig. 

We have explored the parameter space (v, ^') by means of 
an equally spaced grid from v = 3/8tov=lat steps of 1/8, 
and from ^ = 2to*I' = 13 at steps of 0.2. A given model 
will be denoted by the values of the two parameters (v; ^F) in 
parentheses. 

For given (v, T), once the solution for the potential <l)(r) 
is obtained, it can be inserted in the expression of from 
which all the intrinsic and observable profiles and properties 
of the model can be calculated. In particular the density p{r) 
and the anisotropy profile a{r), defined as a(r) - 2 - ((Wg) -H 
require the evaluation of a double integral over the 
velocity space (see also .Bertin & Trenti 2003j . 
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Fig. 1. Relation yiW) for the p^' family of models, for selected 
values of v. To fit the adopted frame, the y values corresponding 
to V = 3/8 have been multiplied by a factor 1/6 and the ones 
corresponding to v = 1/2 by a factor 2/3. 



Fig. 2. Density contrast between the center and the half-mass 
radius along the equilibrium sequence for different values of v. 
The solid line refers to v = 3/8; at fixed ^', models with higher 
values of v are less concentrated. 



3.1. The intrinsic density profile 



The problem of constructing the line profiles F(w,R) for 
a given distribution function of the form f(E, J^) is often 
discussed in the literature (see Gerhard 19911 IGerhard 19931 
ICarol lo et al. 1995 1. The calculation requires the evaluation of 
a triple integral (two integrations over the velocity space or- 
thogonal to the line of sight and one over the radial coordinate 
along the line of sight), which we have performed numerically 
with the same adaptive algorithm used to compute the density 
JBerntsen et al . 19911. In practice, we have followed the proce- 
dure described by Gerhard (1991). The calculation of the quan- 
tities CTproj and hn is then performed as described at the end of 
Sect. 14.31 where they are introduced and defined. 



2.1. The parameter ^ and the density concentration 

In the following we will often refer to the parameter ^' as the 
concentration parameter. Strictly speaking, such term is justi- 
fied only at relatively large values of (*F ^ 5). 

A more intuitive measure of the central concentration of a 
model is given by the ratio piO)/p{rM) of the central density 
to the value of the density attained at the half-mass radius rM- 
As illustrated in Fig. |2l this ratio is a monotonic increasing 
function of ^ only beyond a minimum at^' ^ 4.5. 



3. Density profiles 

In this Section we start with the discussion of the properties of 
the density distribution. 



In the outer parts, i.e. at radii such that r ^ tm, the density 
profile (Fig. |3j is basically the same for all models of the /'"^ 
family. In fact, the dimensionless density associated with the 
f^'^'' distribution function can be expressed, at large values of r, 
as: 



P(r) 



371^ r(2/v) 1 
2 V2 ^ 



(d2(r) + 0(0)3)), 



(5) 



where T{x) is the standard Gamma function 
l |Abramowitz & Stegun 1965) . Curiously, the l/r'* behav- 
ior is found to start from well inside the main body of the 
system, not too far beyond r^, while the density distribution of 
the inner half of the system is similar to that of an isothermal 
sphere. A power-law of the form r fits reasonably well the 
profiles in the transition region from to 3rM. 

As the concentration parameter increases, the point 
where the density profiles merge into a profile common to all 
the models sets in at smaller and smaller radii, so that for in- 
creasing values of T the models appear to converge toward a 
common (singular) model with a central cusp. In these respects, 
the general behavior of the intrinsic density profiles is similar 
to that of the f^o models. Therefore, the behavior of concen- 
trated models is well captured by the following simple formula 
(Uatte 1983> : 



pAr) 



1 



1 



r2 (1 H-r)2' 



(6) 



as shown in Fig. |3] Here, and in the remaining part of this 
subsection, a hat symbol denotes that a quantity is expressed in 
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Log(r/r„) 

Fig. 3. Density profile of some representative /*^' models, with 
V - 3/8, 1/2, and 1, and ^ = 5 and 10. The most concentrated 
model corresponds to (3/8; 10) and the least concentrated to 
(1;5). If the scales are fixed so that = 10 kpc and M = 
lO" Mq, the units for the density are 10"^ M^/pc^. The density 
profiles overlap in the outer parts, beyond a radius that becomes 
smaller and smaller as *P increases. In the plot we also record 
the p J profile, shown here on an arbitrary scale for a convenient 
comparison. 




Fig. 4. Best fit value of the index n (of the R^^" law) associ- 
ated with the projected density profile of the /'^''* models, for 
selected values of v. 



This is in sharp contrast with other modeling procedures that 
are descriptive, rather than predictive (e.g., see lMerritt 1985> . 



a suitable dimensionless form, obviously with no reference to 
the scaling procedure described after Eq. 

In turn, low-*!* models are characterized by a prominent 
core. In fact, some models have density profiles close to that 
of the isochrone model JHenon 19591 1 or of the so-called per- 
fect sphere {pp = (1 + r^)~^). The core may be even represented 
by the density profile of a Plummer sphere iPlummer 19151 1. 

P^'^'^ = (1+^2)5/2 - (7) 

These density profiles are mentioned to better illustrate the 
properties of the mass distribution of the /*^' models in terms 
of well-known profiles. For a discussion of the merits and lim- 
itations of physically based models (such as the /*^' models) 
with respect to other models constructed on the basis of an- 
alytical convenience (e.g., see Merritt 1985 Hernquist 1990 1, 
the reader is referred to the review article by Berlin & Stiavelli 
(1993). On the other hand, in relation to phase-space properties, 
we should emphasize the following point. As will be clear from 
Sect. 4, in the approach adopted in this paper (in which the dis- 
tribution function is constructed from physical arguments) the 
final velocity dispersion and pressure anisotropy profiles can- 
not be set independently, but follow from the self-consistent 
solution. In other words, all the "observable" profiles can be 
seen as consequences of the physical framework considered. 



3.2. The projected density profile 

We have then proceeded to compute a library of projected 
density profiles (which may be compared to observed lumi- 
nosity profiles, under the assumption of a constant M/L ra- 
tio). A first way to characterize these profiles is to fit them 
with the law (Sersic 1968 l. Such fit has been performed 
over a very wide radial range, from O.lRg to IQRg. It shows 
that the /^''' family is well represented by the R^^" law, with 
the index n ranging from 2.5 to 8.5 (see Fig. 0] the slightly 
bumpy behavior of the v = 3/8 curve just reminds us of the 
uncertainties associated with the best-fit determination of n). 
The n - A behavior, characteristic of the de Vaucouleurs law 
(Ide Vaucouleurs 1948> . is mostly associated with concentrated 
high-*? models, but we note that many intermediate-*!' models 
also have the same structural property. The residuals from the 
R^^" best fit (see Fig.|5} are typically within 0.05 mag for con- 
centrated models, while at low values of *!* they are within 0.2 
mag; the general behavior can be compared with that of the /«> 
models (see Fig. A.l in lBertin et al. 2002> . 



4. Phase space properties 

Here we describe the properties of the /'^''* models related to 
the velocity space. 
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Fig. 5. Residuals pi - yu'^", in magnitudes, obtained by fitting 
the law to the projected density profiles of the /*^^ models 
for V = 1 /2 and selected values of *P. The quality of the fit is 
excellent for concentrated models (*P > 7). 



4A. Pressure anisotropy profiles and global anisotropy 
indicators 

The pressure anisotropy of the models can be described by 
means of the anisotropy profile a{r), defined as a(r) - 2 - 
ii^e) + This function, illustrated in Fig.|6j shows 

that the cores are approximately isotropic and that in the outer 
parts the pressure is mostly in the radial direction, in line with 
the qualitative expectations of the violent relaxation scenario 
l |Lynden-Bell 19671 1. Higher values of v are associated with a 
sharper transition from central isotropy to radial anisotropy. 

The degree of anisotropy globally present in our models 
can be characterized in terms of the ratio IKr/Kj, where Kr 
is the total kinetic energy associated with the radial degree of 
freedom, and Kj the corresponding quantity related to the two 
tangential directions. All models present an excess of kinetic 
energy in the radial direction, as illustrated in Fig.0 The ratio 
2K, /Kt is greater than ^1.3 over the whole sequence and be- 
comes larger than 2 in the low-*F region. As a general trend, at 
fixed *F, models with higher values of v are more isotropic. 

In addition to IKrlKj, as a global anisotropy indicator we 
can also refer to the parameter raIrM (see Fig. Q; here the ra- 
dius Fa denotes the anisotropy radius, defined from the relation 
a{ra) = 1. In contrast with the foo models, for which r„ ^ ^ru 
for relatively large values of here we find that concentrated 
models are characterized by ra^rm- 
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Fig. 6. Anisotropy profiles air) of = 5 models for different 
values of v. 



4.2. Velocity dispersion and projected velocity 
dispersion profiles 

As we have seen, the models can be characterized by signifi- 
cant "pressure" anisotropy. This can be illustrated directly by 
the intrinsic velocity dispersion profiles, for which one may 
find significant differences between the tangential and 
the radial dispersion cr^ as far in as r « 0.1 (see Fig.|SJl. 
These velocity space properties, in combination with the den- 
sity distribution, give rise to the projected velocity dispersion 
profiles (calculated from the line profiles, as described in the 
next subsection), which may eventually be compared with the 
observed kinematical profiles; some projected velocity disper- 
sion profiles are shown in Fig.|9l 

4.3. Line profiles 

Pressure anisotropy can affect the shape of the velocity distri- 
bution integrated along a given direction, which can be tested 
observationally by studying the profiles of the lines used to de- 
termine the observed velocity dispersion. In this context, since 
observational limitations prevent us from obtaining accurate 
measurements of line profiles, comparisons between data and 
models are often carried out in terms of certain shape param- 
eters, which measure the deviations from a Gaussian profile 
(e.g., see Gerhard 1991,.de Zeeuw et al. 20 02 1. In view of this 
approach, we have first computed the Une profiles for a number 
of models, following the procedure described in Sect. |5] and 
then we have extracted from them the related value of the h4 
parameter (being non rotating and spherically symmetric, the 
Z'^''* models are associated with line profiles characterized by 
vanishing hi). 
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Fig. 7. Global anisotropy of the /'"^ models shown in terms of 
the parameter IKrjKT (group of curves starting from the top 
left) and of the ratio r^IrM (anisotropy radius to half-mass ra- 
dius; group of curves starting from the bottom left) as a func- 
tion of *F. For low values of *F the global anisotropy IKrjKj 
is basically independent of v. At fixed value of *F the global 
anisotropy decreases with v. 



To extract the velocity dispersion a-pr„j and the value 
hi, we fit the line profile F(w,R), at fixed R, with a 
Gaussian corrected by a fourth order Gauss-Hermite polyno- 
mial ( ,Abramowitz & Stegun 1965) 1: 



F(w,R) = Foexp[-{w/crp,„jY/2]x 

w 



(8) 



{l,,,[12-48(-^)Hl6(-^)^] 

U pro] U pro] 



where Fo, crp,oj, and hi, are free parameters. The fit 
has been performed with a Simulated Annealing method 
(Press et al. 1992 1. 

In general, the deviations from a Gaussian are modest, ^ 
5%, and reduce to within 1% for R ^ R^. The hi parameter 
takes on slightly negative values at the center of the system and 
then becomes positive in the outer parts (see Fig. I10> . in line 
with the results found by Gerhard (1991). 

4.4. The phase space densities N{E, J^) and N{E) 

In view of a comparison with N-body simulations, a physically 
interesting way to characterize the phase space properties of 
the models is in terms of the A^(£', J^) and N{E) densities, 
defined in such a way that M = J N(E)dE = JN(E, J^)dEdfi. 
Therefore, the relation between f^'^XE, J^) and the phase space 
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Fig. 8. Intrinsic "pressure" profiles (cr^/2 - + (w^))/2 

and (T^ = {w^.)) for selected models with v - 1/2. 
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Fig. 9. Projected velocity dispersion profile for selected f'^^^ 
models. If the scales are fixed so that - 10 kpc and 
M - lO" Mq, the velocity dispersion cTproj is given in units 
of 207.4 /fcm/i. 



density N{E, J^) is given by the Jacobian of the transformation 
from d^xd^w to dEdJ^: 



N(E, r) 



2nf^\E, J^) 

Q.r{E,J^) ' 



(9) 
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Fig. 10. Profiles of the parameter for selected /'^^ models 
with V - 3/4 and different values of ^' (the *P = 2 model is 
denoted by squares, *P = 6 by triangles, and *F = 10 by stars). 




Fig. 11. The phase space density N(E, J^) for the (1/2; 6.2) 
model. The graph has been obtained by a Monte Carlo sam- 
pling of the distribution function with 2x10^ points. 



where D.r(E, J^) is the radial frequency of stellar orbits in the 
given potential <l>(r). At variance with the models, the f^'^^ 
models exhibit a singular behavior of N{E, fi) near the origin 
in the dEdJ^ phase space (see Fig. II 1>. 



5. Stability 

The study of the stability of the models could be ap- 
proached either through a linear modal analysis (see Bertin et 
al. 1994) or by means of N-body simulations. Here we prefer 
the latter approach, starting from initial conditions obtained by 



sampling the relevant distribution functions with Monte Carlo 
techniques. 

5.1. General conditions for tine onset of tine radial orbit 
instability 

A simple criterion for the onset of the radial orbit in- 
stability states that instability occurs for systems with 
IKr/Kr ^ 1.75 l|Poiya chenko& ^hukhman 1981) see also 
p-'ridman & Poly achenko 1984 1. We thus note that models with 
*P ^ 4 should be unstable, in a way that is basically indepen- 
dent of V (cf. Fig.Qi. 

5.2. Simulations of stable and unstable configurations 

We have simulated the evolution of different models by 
means of an improved version (see Trenti 2004) of a code writ- 
ten originally by van Albada & van Gorkom (1977; see also 
Ivan A lbada 1982 1. The evolution of the system is collisionless 
because the simulation particles interact with one another via a 
mean field computed with Fourier techniques; the mean density 
associated with the particle distribution is expanded in spher- 
ical harmonics with contributions up to Z = 4 (the code can 
handle contributions up to / = 6). This choice of code is well 
suited for our problem, since its performance has been tested to 
be high for quasi-equilibrium conditions and for smooth spher- 
ically symmetric spatial distributions. To further check our re- 
sults, we have run a comparison simulation with the tree code 
of Londrillo (see ILondrillo et aTT^OOS.) for the (1; 3.2) model 
and found no significant differences. 

In order to sample points from the /'^' distribution func- 
tion and thus to create the initial conditions needed to study the 
stability of the models, we have used an exact inversion 
of the cumulative mass profile and a two-dimensional hit/miss 
method in the (wr,w±) velocity space. The angles needed to 
complete the assignment of the position and velocity vectors of 
the simulation particles are generated randomly. In principle, 
the models extend out to infinity. In practice, we set a cut-off 
radius rcur ~ 50 r^, and we rescale the mass profile in such a 
way that Mivcut) - M. The introduced mass "loss" is generally 
less than 1%. 

Most simulations have been run with 2 x 10^ particles. If 
we refer to a system with total mass M — lO" Mq and half- 
mass radius km - 5 kpc, the evolution is followed from t - 
up to f = te„d = 2 Gyr. This corresponds to more than 20 tj, 
with the natural dynamical time defined as = GM^^^ / {2Kf^^ 
(here K - Kr + Kj is the total kinetic energy); in fact, for the 
range of (v; *!*) considered and the above scales for M and r^, 
we have td ^ (0.65 - 0.80) x lOV- At f ^ lOf^ the system, 
even when it has initially evolved because of initial unstable 
conditions, settles down into an approximate equilibrium state. 

The total energy is typically conserved within 10"^ over 
one dynamical time. The models are non-rotating; the total an- 
gular momentum remains within 10 "* around the small value 
generated in the initialization process. 

To quantify the effects of the radial orbit instability, we 
have especially focused on the evolution of the central den- 
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sity concentration, of the shape of the inertia ellipsoid, com- 
puted on the basis of the distribution of the simulation particles 
within a sphere of radius 3 fm, and of the global anisotropy 
2K,.IKj ratio. The end-products of unstable initial conditions 
are characterized by a final prolate state (with axial ratios 
aa/ai ~ 03/01 < 1), with maximum projected ellipticity con- 
sistent with that of an £3 galaxy. For initially unstable models, 
the central density concentration drops during evolution, but 
the radial scale of the system remains approximately constant. 
For example, for the (1/2; 3) model, the central concentration 
at the end of the simulation is more than 10 times smaller than 
that of the initial conditions, but the radius of the sphere that 
contains 5% of the total mass increases only by 6%, while the 
half-mass radius remains basically unchanged. 

We have studied the evolution of the unstable models in 
terms of exponential growth curves of the form 

^(0 = §o+§i[exp(A:f)-l], (10) 

where go, g\, and the growth rate k are free parameters and g 
represents a global quantity such as the global anisotropy ratio. 
To capture the initial, approximately linear phase of the evolu- 
tion, we have decided to make the fit in terms of the exponential 
growth over a reduced time interval, from t -Qtot - tii„ < t„,^, 
defined implicitly by the relation: 

(f)"'«'4[(f)'-'*(f)'4 

Obviously, this choice of f/,„ is arbitrary and should be consid- 
ered only as a convenient reference time-scale. 

The results of the simulations are reported in Table 1 . The 
threshold for instability is found to be at IKrjKj x 1.70, con- 
sistent with the criterion of Polyachenko & Shukhman (1981). 
Close to conditions of marginal stability, quantities such as 
2K, /Kt vary slowly and by very small amounts, so that the fit 
in terms of exponential growth curves is not well determined. 

An interesting result is the following. We checked whether 
the final state reached as a result of the instability could be 
represented by a model of the family to which the initial 
unstable model belonged. To do this, we fitted the spherically 
averaged density and anisotropy profiles of the end-states of the 
simulations by means of the same family of /^''' models. The 
best fit model thus identified (note the quality of the fit in re- 
lation to both kinematics and density distribution) turns out to 
be the marginally stable model of the sequence with the same 
value of V (see also the interesting arguments and results by 
IPalmeretal. 1990'l l. In the case of v = 1, which we studied in 
greatest detail, the profiles of the end-states obtained starting 
from initially unstable models with ^' in the range from 3 to 4 
are well fitted by the (1;4.2) model (see Fig.ll2>. For the vio- 
lently unstable initial (1; 2) model, the final state is best repro- 
duced by the (7/8; 4) model, which is still moderately unstable. 

We have also run a number of simulations of models ex- 
pected to be stable and checked that indeed the models preserve 
their state for many dynamical time scales. The most concen- 
trated model that we have been able to simulate properly is the 
( 1 /2; 9.4) model, for which we recall that the central concentra- 
tion is p(0)/p(rM) - 1.14 X 10^. This also quantifies the kind of 
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Fig. 12. On the left frames, density and anisotropy profiles 
of the end-products of a simulation starting from the unsta- 
ble (1; 3.2) model (error bars) fitted by the marginally stable 
(1;4.2) model of the /'"^ family (solid line); units for p are 
the units adopted in the simulation. The right frames show the 
residuals Ap = p''*-^* - p and Aa = a*"*-^' - a. 



gradients and ranges of densities that our code is able to handle 
(see Fig.FOt. 



6. A first comparison with tlie observations 

The results in terms of the R^^" law described in Sect. 13.21 
already show that the models possess realistic density 
profiles. This encouraged us to consider a direct compari- 
son with an observed galaxy. For the purpose, we picked the 
round elliptical galaxy NGC 3379, which apparently does not 
possess significant amounts of dark matter (Saglia et al. 1992| 
[Romanowsky et ar^003 note that the models discussed in the 
present paper are one-component models and thus are not ap- 
plicable to systems with prominent dark halos). This galaxy has 
an luminosity profile ( |de Vaucouleurs & Capaccioli 1979| 
[Capaccioli et al. 1990^ . For the kinematics, we considered the 
data of Statler & Smecker-Hane (1999) and recently published 
data-points based on planetary nebulae that extend well beyond 
Re ( [Romanowsky et al. 2003^ . 

The f^^^ model that best describes the data is shown in 
Figs. I14I15I From such model, by adopting a distance of 
1 1 Mpc for the galaxy and an absolute magnitude in B band 
Mb = -20.0, we obtain a mass-to-light ratio M/Lb = 4.7 in 
solar units. Population synthesis models for NGC 3379 predict 
a mass-to-light ratio between 4 and 9 (IGerhard et al. 200ll . In 
comparison, Romanowsky et al. (2003) report a mass-to-light 
ratio M/Lb = 7.1 + 0.6. 
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Fig. 13. Evolution of the fractional mass radii (of the sphere 
containing, from the bottom to the top, 6.25 x 10"^, 1 .25 x 10""*, 
0.01, 0.05, 0.5, 0.7, 0.95, 0.99 M, respectively) for a simulation 
with 8x10^ particles of the concentrated model (1 /2; 9.4). 
The model is stable over the time interval considered (about 10 
dynamical times). The radius r and the time t are given in units 




of 10 kpc and 10** yr. 



l.og(tx/Rj 

Fig. 14. Comparison between the photometric profile of NGC 
3379 and the (l/2;9.4) model of the family. The photo- 
metric profile in B is taken from de Vaucouleurs & Capaccioli 
(1979). Other members of the family (e.g. the (1;9.2) model) 
could better reproduce the observed data especially at large 
radii, but they would perform less well on the velocity disper- 
sion profile. 



Table 1. Growth rate k, with an estimate of the related uncer- 
tainty /Sk, for the radial orbit instability from the simulation of 
a set of models. Here k - IKrlKj is the global anisotropy 
parameter and rj - 03/01 is the axis ratio of the inertia ellipsoid 
referred to the sphere of radius 3 rM- For M - 10" M© and 
= 5 kpc, k is given in units of 10"^ yr^^ . 



(v;1') 


k 


AA: 




K(te„d) 




(1;2.0) 


0.85 


0.1 


2.65 


1.87 


0.65 


(1;3.0) 


0.48 


0.07 


2.17 


1.84 


0.67 


(1;3.2) 


0.37 


0.07 


2.09 


1.84 


0.69 


(1;3.4) 


0.34 


0.07 


2.02 


1.83 


0.74 


(l;3.6) 


0.30 


0.05 


1.95 


1.85 


0.80 


(l;3.8) 


0.16 


0.04 


1.89 


1.82 


0.85 


(l;4.0) 


0.08 


0.05 


1.84 


1.81 


0.91 


(1;4.2) 


0.007 


0.10 


1.78 


1.76 


0.94 


(l;5.0) 


0.001 


0.01 


1.60 


1.60 


0.99 


(l;9.0) 


< 10" 




1.32 


1.32 


0.99 


(3/8;3.0) 


0.90 


0.10 


2.28 


1.86 


0.73 


(3/8;5.0) 


0.45 


0.15 


1.75 


1.71 


0.93 


(l/2;3.0) 


0.70 


0.10 


2.21 


1.86 


0.70 


(l/2;4.0) 


0.43 


0.05 


1.92 


1.80 


0.85 


(l/2;5.0) 


0.1 


0.2 


1.68 


1.67 


0.96 


(l/2;6.0) 


< lo - '' 




1.54 


1.54 


0.99 
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0.64 
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1.85 


0.70 
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0.16 


0.02 


1.86 
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< 2 10^^ 




1.63 


1.63 
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Fig. 15. Kinematic data for NGC 3379 described in terms of 
the (l/2;9.4) model of the family. The inner data (plain 
eiTor bars) are taken from the stellar spectroscopy of Statler & 
Smecker-Hane (1999), while the outer data (crosses with error 
bars) refer to the binned velocity dispersion determined from 
the study of planetary nebulae ( |Romanowsky et al. 2003| l. 
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7. Discussion and conclusions 

In this paper we have studied the structural and dynamical 
properties of a two-parameter family of models of partially re- 
laxed stellar systems. These models had been proposed earlier 
JStiavelli & Berlin 1987 1 as physically relevant to the galaxy 
formation scenario based on collisionless collapse and incom- 
plete violent relaxation. They had been noted to possess some 
realistic properties (for v in the range 1/2-1 and for relatively 
large values of *!*). However, they had been left basically aside 
and not studied further in systematic detail. Additional physical 
interest was noted recently l lBertin & Trenti 200 3 ), in an inves- 
tigation focused on their thermodynamic properties, in relation 
to the paradigm of the gravothermal catastrophe. Because of 
such physical interest, we have decided to undertake a thor- 
ough comparison between the models and the products of col- 
lisionless collapse, as generated in N-body simulations. Such a 
comparison will be the subject of a follow-up paper, of which 
the present article forms the necessary basis. In that paper we 
will consider a relatively large set of simulations of collision- 
less collapse, starting from a variety of initial conditions. We 
will address the issue of the conservation of Q and identify 
the range of v for which the global quantity is best conserved. 
Furthermore, we will show that in many cases the f-^^ mod- 
els can provide a surprisingly good fit to both the density and 
the pressure profiles, over nine orders of magnitude of the den- 
sity distribution (see lBertin & Trenti 2004t . The best-fit models 
will turn out to be close to marginal stability with respect to the 
radial orbit instability. 

Here we have shown that the family of /^''' models ex- 
hibits a variety of structural properties, within a common gen- 
eral behavior. The model characteristics can be summarized by 
referring to three separate regimes: low-*? models (typically, 
*P < 4), intermediate-*!' models (typically, 4 < < 8), and 
high-*P models (typically, > 8). In practice, the values of 
that mark the transition between different regimes do not de- 
pend significantly on the value of v, at least in the range ex- 
plored in this paper (i.e., from v = 3/8 to v = 1). 

In the intermediate-*!' and the high-*!' regimes the concen- 
tration ratio p(Q)lp(rM) increases monotonically with *!'; mod- 
els with lower values of v have larger values of p{0)lp{rM) (by 
up to one order of magnitude at given *!*, in the explored range 
of v). In the same regimes, as *!' increases, the density distri- 
bution converges onto a common profile characterized by an 
approximate r^'^ behavior at r > and by an approximate r^^ 
behavior at r < r^, with an inner core that becomes smaller 
and smaller. In contrast, models in the low-*!' regime have a 
substantial core, comparable in structure to that of a Plummer 
or of an isochrone model. In general, the projected density dis- 
tribution of the models is well fitted by the law, with n 
ranging from 2.5 to 8.5. Typically, high-*!' models are accu- 
rately described in terms of the R^^"^ law. 

The models are all characterized by significant radially- 
biased pressure anisotropy. In terms of the global anisotropy 
ratio IKr/Kr, at given *!' models in the low-*!' regime have sim- 
ilar amounts of pressure anisotropy, which increases rapidly as 
*!* decreases. In practice, all low-*? models should all be unsta- 
ble with respect to the radial-orbit instability. Intermediate-*!' 
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and high-*!' models appear to be safely stable, but they may 
be only marginally so for the lowest value of v that we con- 
sidered. This latter remark also explains why we have decided 
not to move below v = 3/8, given the fact that radial pressure 
anisotropy is larger for lower-v models. In terms of velocity dis- 
persion profiles, to some extent pressure anisotropy is already 
significant at r 0.1 r^, even for stable models. However, out 
to r X tm the line profiles of stable models deviate very little 
from a Gaussian, with /14 within the 5% level. 

By means of an improved version of a code introduced 
earlier by van Albada, we have checked the properties of the 
models with respect to the radial-orbit instability, confirming 
the general expectations based on the global anisotropy pa- 
rameter 2K,-/Kt- Curiously, we find that unstable models tend 
to evolve by staying close to the equilibrium sequence, while 
moving up in *!' so as to reach marginal stability. Furthermore, 
we have been able to show that the N-body code is able to 
handle highly concentrated models (up to *!' = 9.4, with 
p(Q)lp(rM) - 1.14 X 10^), of which it has demonstrated the 
long-term stability. 

Finally, a comparison with the observed surface brightness 
and kinematical profiles for the galaxy NGC 3379 has shown 
some merits and some limitations of the one-component fam- 
ily of /*^* models when applied to observed objects. In partic- 
ular, it is interesting to see that, within the extremely idealized 
framework at the basis of the construction of these models (see 
also the general comments made in Sect. 1.1), the observed lu- 
minosity profile is very well reproduced over about ten magni- 
tudes; at the same time the models can well reproduce the inner 
part of the relevant kinematical profile, inside . This compar- 
ison supports the view that the models studied in this paper 
may be helpful to describe the luminous component of elUp- 
tical galaxies. The obvious limitations of the one-component 
models that we are considering are clearly brought out by their 
inadequacy to capture the change in the observed kinematical 
profile occurring around and thus their failure to reproduce 
such profile in the outer parts. This is interpreted as a signature 
of the presence of an additional component that our models, in 
the form developed so far, a priori ignore. 
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Londrillo, for their collaboration and for a number of useful comments 
and suggestions, and N. Douglas and A. Romanowsky, for providing 
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contact M. Trenti, who will be glad to provide the relevant numerical 
routines to initialize N-body simulations and to compute intrinsic and 
projected profiles for the family of models. 

Appendix A: Constants and scales 

From the definition of the /*^*'' distribution function (see Eq.|2j, 
the constants (A, a, d) have the following dimensions: 

A = [ML-^T\ (A.l) 
a = [L-^T\ (A.2) 
d = [L-^'^T-'''^ (A.3) 

Therefore, as scales for mass and radius we may refer to: 



12 



M. Trenti & G. Berlin: A family of models of partially relaxed stellar systems 



(A.4) 
(A.5) 



Then we introduce the dimensionless quantities r = r/Rscaie 
(and thus = /m(v,^) = ru/Rscaie), M = M(v^^') = 
M/Mscaie, where M(v,*P) = j pir)d^f, w = o'^^vv, and <J) = 
for radius, mass, velocity, and potential, respectively. 

Recalling the expression for y (see Sect. 2, after Eq. (|2jl), 
for given v we can express (A, a, t/) in terms of (M, r^, *I0- 



1 



Tm M(y,>P) 



47rGr(v, >F) rM(v, "¥) M 



d = [4■nGy(v,^')] 



A = 



v/4 



pM(v,'P) 



\M(V,>1')/ \ f-M 



5 v/4 



[47rGr(v,>F)]3/2 \ 



rM(v,*l')\^^VM(v,>F)^'^^ 
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